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The classic metallurgical systems - noble metal alloys - that have formed the benchmark for 
various alloy theories, are revisited. First-principles fully relaxed general potential LAPW total 
energies of a few ordered structures are used as input to a mixed-space cluster expansion calculation 
to study the phase stability, thermodynamic properties and bond lengths in Cu-Au, Ag-Au, Cu- 
Ag and Ni-Au alloys, (i) Our theoretical calculations correctly reproduce the tendencies of Ag-Au 
and Cu-Au to form compounds and Ni-Au and Cu-Ag to phase separate at T — K. (ii) Of all 
possible structures, CU3AU (LI2) and CuAu (Llo) are found to be the most stable low-temperature 
phases of Cui-^ Au^ with transition temperatures of 530 K and 660 K, respectively, compared to the 
experimental values 663 K and w 670 K. The significant improvement over previous first-principles 
studies is attributed to the more accurate treatment of atomic relaxations in the present work, (iii) 
LAPW formation enthalpies demonstrate that LI2, the commonly assumed stable phase of CuAu3, 
is not the ground state for Au-rich alloys, but rather that ordered (100) superlattices are stabilized, 
(iv) We extract the non-configurational (e.g., vibrational) entropies of formation and obtain large 
values for the size mismatched systems: 0.48 fcs/atom in Nio.sAuo.s (T = 1100 K), 0.37 fcs/atom in 
Cuo.141Ago.859 (T = 1052 K), and 0.16 fc s /atom in Cuo.sAuq.s (T = 800 K). (v) Using 8 atom/cell 
special quasirandom structures we study the bond lengths in disordered Cu-Au and Ni-Au alloys 
and obtain good qualitative agreement with recent EXAFS measurements. 

PACS numbers: 61.66.Dk, 71.20.Gj, 81.30. Bx 



I. INTRODUCTION: CHEMICAL TRENDS IN 
NOBLE METAL ALLOYS 

Noble metal alloys are, experimentally, among the 
most studied intermetallic systems. 1-24 In addition, 
the Cu-Au system has been considered the clas- 
sic paradigm system for applying different theoret- 
ical techniques of phase diagram and phase stabil- 
ity calculations. 25-63 Most notably, this system has 
been considered as the basic test case for the clas- 
sic Ising-hamiltonian statistical-mechanics treatment of 
alloys. 25-32 More recently, noble metal binary alloys have 
been treated theoretically via empirical fitting of the 
constants in Ising hamiltonians, 25-34 semiempirical in- 
teratomic potentials, 35-47 and via first-principles clus- 
ter expansions. 48-55 The essential difference in philoso- 
phy between the classical application of Ising models to 
CuAu 25-30 ' 33 and more modern approaches based on the 
density functional formalism 6 is that in the former ap- 
proach the range and magnitudes of the interactions are 
postulated at the outset (e.g., first or second neighbor 
pair interactions), while the latter approaches make an 
effort to determine the interactions from an electronic 
structure theory. However, despite recent attempts, 48-54 
it is still not clear whether the noble metal alloys can 
be essentially characterized as systems with short-range 
pair interactions, or not. 

Now that first-principles cluster expansion 
approaches 65 ' 66 have advanced to the stage where both 



T = ground state structures and finite-temperature 
thermodynamic quantities can be predicted without any 
empirical information, it is interesting to take a global 
look at the noble metal alloy family. Table I summarizes 
some of the salient features 1-4,14 ' 15 ' 18 ' 67-69 of the four 
binary systems Cu-Au, Ag-Au, Cu-Ag and Ni-Au. We 
included the relative lattice constant mismatch Aa/a = 
2 \cla — clb\ I \ a A + Ob| between the consituents, 67 the 
electronegativity difference A% = xa~Xb on the Pauling 
scale, 68 the mixing enthalpy of the equiatomic alloy, 2,18 
the sign of the calculated nearest neighbor pair interac- 
tion J2 (present study), the structural identity of the low- 
temperature phases 1-4,67 and the order-disorder transi- 
tion (or miscibility gap) temperatures 2,69 T c . Some in- 
teresting observations and trends which we will attempt 
to reproduce theoretically, are apparent from this general 
survey: 

(i) Despite a large (12%) size mismatch in Cu-Au, and 
a small (~ 0%) size mismatch in Ag-Au, both systems 
form ordered compounds at low temperatures and have 
negative mixing enthalpies, suggesting attractive ("anti- 
ferromagnetic" ) A-B interactions. Thus, when the differ- 
ence in the electronegativity A% of the constituent atoms 
is sufficiently large, as it is in CuAu and AgAu, size mis- 
match apparently does not determine ordering vs. phase 
separation tendencies. 

(ii) Despite a similar size mismatch (12%) in Cu-Au 
and Cu-Ag, the former orders while the latter phase- 
separates. Thus, the existence of large electronegativity 
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TABLE I. Major physical properties of Ag-Au, Cu-Ag, Cu-Au and Ni-Au alloys. We give constituent size mismatches, 
Aa/a — 2(aA — aB)/(aA + fflB), electronegativity differences on the Pauling scale, 68 Ax, mixing enthalpies of the disordered alloys 
at the equiatomic composition, A_ff m i x (a; = signs of the nearest-neighbor Ising interaction, J2, order-disorder transition 
temperatures (or miscibility gap temperatures for Cu-Ag and Ni-Au), T c (x — |), and excess entropies of solid solutions, 
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difference in Cu-Au (as opposed to the small difference 
in Cu-Ag), seems to induce ordering tendencies. 

(iii) Cu-Ag and Ni-Au both phase-separate (and have 
positive Ai? m ; x ) as they have large size mismatches. 
Yet, Ni-Au having a large electronegativity difference, 
shows an ordering-type nearest-neighbor pair interaction 
( J2 > 0), just like the compound forming Cu-Au and Ag- 
Au, while Cu-Ag has a clustering-type nearest-neighbor 
interaction (J 2 < 0). Thus, the sign of J 2 does not reflect 
the low temperature ordering vs. phase separation. 

(iv) The amount ASxs = AS^* - A5 idc ai 
by which the measured entropy 2 A5t* t pt deviates 
from the ideal configurational entropy ASidcai = 
fee log a; + (1 — x) log(l — x)], is unexpectedly large in 
Cu-Ag and Ni-Au, indicating a large non-configurational 
entropy of formation. 

Other interesting facts about the noble metal binary 
intermctallics include: 

(v) Despite numerous studies, 1-4 ' 7 ' 8 ' 10 " 12 the struc- 
ture of the ordered phases in Au-rich Cu-Au is not well 
established yet. It is often assumed 1-4 that the sta- 
ble Au-rich low-temperature phase is CuAu3 in the Ll 2 
structure, but direct experiments 7 ' 8 ' 10 below the order- 
disorder transition temperature T c (x = |) w 500 K 
are difficult because the diffusion rates are very low and 
even the best ordered samples contain significant disor- 
der. Possible further thermodynamic transformations at 
lower temperatures may be kinetically inhibited. 

(vi) The trends in bond lengths vs. composition 
are non-trivial. Traditionally, all coherent-potential- 
approximation based theories 70-72 of intermetallic alloys 
have assumed that the nearest-neighbor bond lengths are 
equal, Raa = Rab = Rbb, an d proportional to the av- 
erage lattice constant. Recent theories 73-75 suggested, 
however, that bond lengths relax in the alloy to new val- 
ues, and this has a significant effect on the electronic 
structure. 53 ' 76 ' 77 Recent EXAFS experiments on NiAu 23 



and CuAu 24 show distinct Raa 7^ Rab 7^ Rbb bond 
lengths, which need to be explained. 

In this work we will analyze the above mentioned 
trends in terms of a first-principles mixed-space cluster 
expansion, 65 ' 66 based on modern local density approxi- 
mation (LDA) total energy calculations. We reproduce 
the observed trends (i)-(vi) in ordering preferences, mix- 
ing enthalpies A_ff m j x , transition temperatures T c and 
interatomic bond lengths. In addition, we predict new, 
hitherto unsuspected ordered phases in Au-rich Cu-Au 
alloys. 

II. BASIC IDEOLOGY AND METHODOLOGY 

There are many problems in solid state physics that 
require knowledge of the total energy E{a) of a lattice 
with N sites as a function of the occupation pattern a 
of these sites by atoms of types A and B. This informa- 
tion is needed, for example, in the ground state search 
problem, 72 where one seeks the configuration with the 
lowest energy at T = K. {E(a)} is also needed for calcu- 
lating the temperature- and composition-dependent ther- 
modynamic functions and phase diagrams of an A\_ X B X 
alloy. 

A direct, quantum- mechanical calculation of the to- 
tal energy E dilcct (a) = (&\ H\^}/(^\^} (where * is the 
electronic ground state wave function and H is the many- 
electron Hamiltonian) is possible only for a limited set 
of configurations a. This is so because (i) the com- 
putational effort to solve the Schrodinger equation for 
a single configuration scales as the cube of the number 
of atoms per unit cell, so that only small unit cells can 
be considered, (ii) there are 2 N configurations, and (iii) 
for each configuration, one has to find the atomic relax- 
ations Su min (a) which minimize the total energy. Conse- 
quently, one searches for a "cluster expansion" (CE) that 
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accurately reproduces the results of a direct, quantum- 
mechanical (e.g., LDA) calculation 

E CE (<j) = £direct(<7), (1) 

without the unfavorable scaling of the computational cost 
with the size of the unit cell. 

In designing a cluster expansion, there are few choices 
of independent parameters. For example, one could 
choose to obtain a cluster expansion for the volume- 
(V) dependent equation of state indirect V) [see, e.g., 
Refs. 52, 78, 79], or to find a cluster expansion for the en- 
ergy at the volume V min (a) that minimizes indirect (c, V). 
We choose the latter possibility. Furthermore, for each 
configuration a, we wish to reproduce the total energy 
corresponding to the fully relaxed cell shape and atomic 
positions {<5u m i n (cr)}. In other words, we choose to rep- 
resent 

Ece(v) = E dilcct [cr; Su mm (a);V mlD (a)} = E dilcct (a). (2) 

Note that by focusing on the equilibrium energy of each 
configuration, we give up the possibility of studying non- 
equilibrium geometries (e.g., bond lengths) and equations 
of state. Instead, for each occupation pattern a, we can 
find the total energy E(a) of the atomically relaxed and 
volume-optimized geometry. 

The best-known cluster expansion is the generalized 
Ising model in which the equilibrium total energy of an 
arbitrary configuration a is expanded in a series of basis 
functions defined as pseudospin products on the crystal 
sites: 

E(a) = Jq + JiSi + 2 X/ Jij^iSj 

i ijtj 
+ y ' JijkSiSjSk + . . . , (3) 

where in binary A\- X B X alloys Si = +1 or —1, depending 
on whether the site i is occupied by an atom of type A 
or B. Equation (3) is valid whether the lattice is relaxed 
ot not, as long as a one-to-one correspondence exists be- 
tween the actual positions of atoms and the ideal fee 
sites. The practical usefulness of the cluster expansion 
Eq. (3) rests on the assumption that the effective cluster 
interactions (ECI's), J^ , Jijk, ■ ■ ■ , are rapidly decreasing 
functions of the number of sites and intersite separation, 
so that only a finite number of terms has to be kept in 
Eq. (3) for the desired accuracy. In this case, we can 
write the formation enthalpy of structure a, 

AH diiect {a) - E{a) - xE A - (1 - x)E B , (4) 

where Ea and Eb are total energies of the pure con- 
stituents A and B, as the following cluster expansion 
(CE): 

N f 

AH CE (a) = Jo + J2 D f J f U f(°)- ( 5 ) 

/ 



Here Nf is the number of nonzero effective interactions 
and Hf(a) are lattice averages of the spin products in 
configuration a. 

Sanchez, Ducastelle and Gratias 80 have shown that 
there is a set of composition-independent interactions 
for Eq. (3) which can exactly reproduce the directly 
calculated total energies of all configurations a. This 
statement is strictly true if all possible clusters are in- 
cluded in Eq. (3), and should hold for the truncated se- 
ries Eq. (3) if the cluster expansion is well converged. 
Several methods 81,82 yield concentration-dependent ef- 
fective interactions, providing in principle equally valid 
schemes for representing AiJdiroct (c) in terms of a cluster 
expansion. In the present work, we select composition- 
independent interactions, since these can be directly com- 
pared to previous Ising model treatments 25 " 34,48-55 of 
the noble metal alloy phase diagrams. 

A number of issues arise in trying to construct a cluster 
expansion that satisfies Eq. (2): 

(i) The number of interactions and their types (pair, 
multibody) cannot be decided arbitrarily, but must be 
constrained by a microscopic electronic-structure theory 
according to Eqs. (1) and (2). 

(ii) In most configurations a, atoms move away from 
the ideal lattice sites, which not only lowers the total en- 
ergies -Edirect(c), but also slows down the convergence 85 
of the expansion Eq. (3). The solution is to have a clus- 
ter expansion with many interaction terms Nj that can 
represent such situations. We accomplish this by using 
a reciprocal space formulation, formally equivalent to an 
infinite number of real-space pair interactions. 

(iii) Some cluster expansions 78 require that the number 
of interactions Nj must equal the number of configura- 
tions N a whose total energies need to be evaluated via 
the direct electronic-structure method. The number of 
such calculations may be excessive in view of (ii). We 
thus introduce a method in which N a <C Nj. Further- 
more, interactions that are not needed to satisfy Eq. (2) 
are automatically discarded. 

(iv) One has to deal with the macroscopic elastic strain 
leading to a k — > singularity in the Fourier transform 
of the pair interactions, 

J pair (k) = ^V P air(R, - R,)e- jkR ^ , (6) 

j 

where Dj is the number of {Rj, Rj} pairs per lattice site. 
As shown by Laks et al. 65 (see also the discussion below), 
in size mismatched systems the correct J pa ; r (k) depends 
on direction k in the long-wavelength limit k — > 0. To 
solve this, we express J pa i r (k) as a sum of two parts, 

Jpair(k) = J SR (k)+ J CS (fc), (7) 

where </sR,(k) is an analytic function of k and can be 
obtained from short-ranged real space pair interactions, 
while Jcs(k) contains the nonanalytic behavior around 
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k = and depends only on the direction k. To ex- 
plain this singularity, we consider the energy of a co- 
herent A n B n superlattice, formed by a periodic stacking 
of n layers of A and n layers of B in direction G. By 
introducing the structure factor, 



(8) 



the total pair interaction energy in Eq. (3) can be ex- 
pressed as a reciprocal space sum: 



£ pair (<7) =^Jpair(k) \S(k,a)f 



(9) 



A n B n superlattice has nonzero structure factor at 
k = j-G, and therefore its energy is determined by 
Jpaii(^G). As n — ► oo its formation energy is given 
by a sum of the epitaxial deformation energies of pure 
constituents needed to bring them to a common lattice 
constant in the plane perpendicular to G. Since the epi- 
taxial deformation energy of pure constituents is direc- 
tion dependent (e.g., it is easier to stretch Cu in [100] 
planes than in [111] planes, see Sec. IIIB), the formation 
energy /S.H{A 00 B 00 ) is also direction dependent. There- 
fore, linik^o Jp a i r (k) must depend on the direction of 
approach to the origin, proving that J pa i r (k) is singu- 
lar. Physically, the nonanaliticity of J pa i r (k) is caused 
by long-range interactions via macroscopic elastic strain 
and cannot be reproduced using finite-ranged real-space 
pair interactions, but must be accounted for explicitly in 
reciprocal space. If the singularity is neglected, then as 
explained in Ref. 65, the cluster expansion fails not only 
for long-period (n — > oo) superlattices A n B ni but also 
for those short-period (n > 2) superlattices which have 
not been explicitly included in the constraint Eq. (2). We 
emphasize that although the contribution of J s ; ng (k) to 
the formation energy is nonzero only in size-mismatched 
systems, it is not related to the atomic relaxation energy 
for a particular structure a in any simple way (except if 
a itself is a long-period superlattice). 

The singularity in J pa j r (k) is similar to the singularity 
in the dynamical matrix D a p(KK'\k) of polar crystals in 
the long-wavelength limit, 83 caused by long-range elec- 
trostatic interactions via macroscopic electric field. In 
lattice dynamics, D a p (kk' |k) is expressed as a sum of 
regular and singular parts, D a p (kk! |k) = D s ^ g (KK'\k) + 
D T °J*(KK'\k), where D^f(/t/t'|k) (analytic as k -> 0) is 
due to short-range force constants. The singular part 
D™ s (KK'|k) gives rise to LO/TO splitting of the zone- 
center optical frequencies lot in polar crystals, and also 
leads to a directional dependence of u>r(k) in uniaxial 
crystals (e.g., CuPt-type GaInP2). These phenomena 
cannot be reproduced by any set of finite-ranged micro- 
scopic force constants, but have to be calculated explic- 
itly using the macroscopic Maxwell equations. 84 

In summary, we seek to find a function _Ece(c) 
which accurately reproduces the LDA total energies 



E LDA [a, <5u min (cr); V m i n (a)} = E hDA (a) at the atomi- 
cally relaxed geometry and equilibrium volume of con- 
figuration a. The function Ece(p) we consider includes 
composition- and volume-independent interactions, so as 
to maintain maximum similarity with the classical Ising 
model. The number and type of interactions are not de- 
cided arbitrarily, but are constrained by the electronic 
structure theory used (here, the LDA). Relaxation is 
treated accurately by including long-range pair interac- 
tions in the reciprocal space representation. The k — > 
singularity, affecting both short and long-period super- 
lattices, is dealt with explicitly. 

The above requirements are satisfied by the mixed 
space cluster expansion (MSCE): 



Aff C E(a)=^J pair (k) |S(k,a)| s 



MB 



J2DfJfU f (a) + AEcs(cT). 
f 



(10) 



We have separated out the so-called equilibrium con- 
stituent strain energy term, AEcs(°~), which accounts 
for the k — > singularity. 65 In Eq. (10) we do not need 
to calculate AEqs(o~) for each configuration a, but only 
for the directions k of the wave vectors with nonzero 
S(k, a). In fact, it is constructed to coincide with the 
elastic strain energy of coherent superlattices in the long- 
period limit: 65 



A£; cs (a)-^J cs (x,fc)|5(k,a)| 2 , 

k 



(ii) 

(12) 



where S(k, a) is the structure factor from Eq. (8). The 
quantity AE^? s (x,k) depends only on the direction k, 
and will be given in Sec. IIIB. Equation (11) is ex- 
act for long-period superlattices, but represents a choice 
for short-period superlattices and non-superlattice (e.g., 
Ll 2 ) structures. It has been found 65 that the choice 
Eq. (11) improves the cluster expansion predictions also 
for short-period superlattices. 

Equation (10) is a generalized Ising model description 
of the formation energy of any relaxed configuration a, 
even if a direct LDA calculation for this a is impractical. 
The cluster interaction energies {J pa ir(k)} and {J/} are 
obtained by fitting Eq. (10) to the LDA formation ener- 
gies. An additional smoothness requirement is imposed 
on J pa ; r (k), which ensures that the pair interactions are 
optimally short-ranged in real space. Namely, we mini- 
mize the sum 

AL = ^~ E w ° l AH Cv(<?) - AH LDA (a)f 



+ ^ E J P-W [- V k] V2 J P-(k), (13) 
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TABLE II. Definition of the small-unit-cell ordered structures used in the LDA total energy calculations. 
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where A and t are free parameters and a is a normaliza- 
tion constant. 65 Typically we choose A = 4 and t = 1, 
but the fit is not sensitive to this choice. 

This approach solves the four problems indicated above 
in the sense that (i) the fitting process itself automati- 
cally selects the pair interactions that are essential to 
obtain a good fit (process still does not select multibody 
figures), (ii) the pair interactions can be of arbitrary long 
range, facilitating treatment of systems with large elastic 
relaxations, (iii) the number of pairs can be much larger 
than the number of ordered structures in the fit, and (iv) 
the directly calculated constituent strain energy AEcs 
contains the k — > singularity. Unlike all CPA-based 
methods, 70,71 the present approach includes full account 
of atomic relaxation and local environment effects. Un- 
like the classical Ising descriptions, 25 ' 27-33 the interac- 
tion energies are determined by the electronic structure 
rather than being guessed. Finally, unlike the compu- 
tational alchemy linear response approach, 85 multibody 
terms are included here. 

Having written the expression for the total energy of 
arbitrary configuration, Eq. (10), we can evaluate its con- 
stants from a limited number of LDA calculations on 
small unit cell (A^ atom s < 10) ordered structures with 
fully relaxed atomic positions. Equation (10) can then 
be employed in simulated annealing and Monte Carlo 
calculations 86 ' 87 yielding T — ground states and T > 
statistical and thermodynamic properties. Further de- 
tails of the method are given in Sec. III. 



III. DETAILS OF THE METHOD 

A. T = energetics 

The calculations of T = total energies employ the 
full-potential linearized augmented plane wave method 88 
(FLAPW). The basis set consists of plane waves in the 
interstitial region, augmented in a continuous and differ- 
entiable way with the solutions of the radial Schrodingcr 
equation inside the non-overlapping muffin-tin spheres. 
Non-spherical potential and electronic charge density 
terms are calculated in all space and included in the 
Hamiltonian matrix. Core states are treated fully rela- 
tivistically and recalculated in each self-consistency iter- 
ation. The wave equation for the valence states includes 
all relativistic effects except the spin-orbit interaction, 
i.e., they are treated scalar relativistically. FLAPW is 
the most accurate all-electron method, superior to the 
methods employing overlapping atomic spheres (atomic- 
spheres approximation - ASA) and/or shape approxima- 
tions to the potential. 

We use the Wigner exchange-correlation functional. 89 
As a check, we have performed several calculations using 
the Perdew-Zunger 90 parametrization of the Ceperley- 
Alder 91 functional and the generalized gradient approxi- 
mation of Pcrdew and Wang. 92 We find (see Sec. IV A 1) 
that the various exchange-correlation functional change 
the enthalpies of formation of ordered Cu-Au compounds 
by a negligible amount (less than 2 meV/atom). 
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TABLE III. LDA calculated formation [Eq. (4)] enthalpies for fee superstructures (defined in Table II) of Ag-Au, Cu-Ag, 
Cu-Au and Ni-Au. The numbers in parentheses represent errors of the cluster expansion fit. All energies in meV/atom. 
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The total energies of the ordered structures and end- 
point constituents are obtained keeping all computa- 
tional parameters exactly equal. Specifically, we always 
use the same basis sets (RK mayi = 9), charge density 
cutoffs (iJ.K'max = 19), muffin-tin radii Rau = 2.4ao, 
i?Ag = -Rcu = Rm = 2.2ao, maximum difference in the 
angular momenta in the nonspherical Hamiltonian terms 
(^max = 4) , maximum angular momenta in the nonspher- 
ical charge densities and potentials inside the muffin-tin 
spheres (7 max = 8), and equivalent k point sets 93 in the 
evaluation of Brillouin zone integrals. When the unit cell 
vectors of the ordered compound permit, we choose a k 
mesh equivalent to the 60 special points 8 x 8 x 8 fee 
mesh. Several structures (e.g., those of A 2 B or AB 2 sto- 
ichiometry) have reciprocal unit cell vectors which are 
incommensurate with the 8x8x8 mesh. In these cases 
we calculate the total energies of the compounds and fee 
constituents with a finer k point grid. This procedure 
ensures that, due to systematic cancellation of errors, 
the formation enthalpies AH(a), Eq. (4), converge much 
faster than the total energies. Indeed, the tests for Cu- 
Au described in Sec. IV A 1 show that with our choice of 
parameters AH(a) are converged to within 2 meV/atom. 

The atomic positions are relaxed using quantum 
mechanical forces 94 obtained at the end of the self- 
consistency iterations. Minimization of the total energy 
with respect to the cell-external degrees of freedom is 
done by distorting the shape of the unit cell and tracing 
the decrease in the total energy. We estimate that the for- 
mation enthalpies are converged to at least 5 meV/atom 
with respect to all relaxational degrees of freedom. 

Table II and its caption defines our small-unit-cell or- 
dered structures. Many are actually superlattices along 
(100), (110), (111), (201) and (311) directions. Table III 
gives the calculated LDA formation energies [Eq. (4)] for 
these Au-Ag, Cu-Au, Cu-Ag and Ni-Au compounds. 
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FIG. 1. A schematic illustration of the concept of the epi- 
taxial softening function q(G), given by the ratio of the bulk 
(upper curve) and epitaxial (lower curve) deformation ener- 
gies. In the harmonic approximation q(G) is the ratio of the 
curvatures of these curves at the equilibrium point. 



energies required to deform bulk A and B to the epitaxial 
geometry with a planar lattice constant a s : 



AE° c q s {x, G) = min xAE A p \a s , G) 

Q-s L 

+ (l-x)AE^(a s ,G) 



(14) 



Here A_E cpl (a s , G) is the strain energy of the material 
epitaxially stretched to the lattice constant a s in the di- 
rection orthogonal to G, and then allowed to relax along 
G. AE cpl (a s , G) is related to the bulk equation of state 
Ai? bulk (a s ) via the epitaxial softening function q(a s ,G): 



B. The constituent strain energy 



q(a s ,G) 



AE^jas^) 
A£ bulk (a s ) : 



(15) 



It is well known 66 that real-space cluster expansions 
with finite-ranged interactions incorrectly predict zero 
formation enthalpies per atom for coherent long-period 



where Ai?^ (a s ) is the energy required to hydrostati- 
cally deform a bulk solid to the lattice constant a s . Fig- 
ure 1 illustrates the concept of epitaxial softening: 97 when 



A p B q superlattices, while the correct answers are non- the bulk solid is deformed hydrostatically from a cq to 

a s a eq , its energy rises. Energy can then be lowered if 
we keep a x = a y = a s but relax the third lattice vector 
to its equilibrium value. q(a s ,G) measures the relative 
energy lowering. 

Figure 2 shows the calculated LDA q's for Cu, obtained 
by minimizing the total energy with respect to the lat- 
tice constant c parallel to G for each value of the sub- 
strate lattice parameter a s . As explained in Ref. 96, this 
treatment neglects the so-called shear strain terms, but is 
exact for the high symmetry directions (100), (111) and 
(110). The calculated qcu(a s , G) is seen to be a nontrivial 
function of the substrate lattice parameter a s and direc- 
tion G. In contrast, the harmonic elasticity theory, 97-100 



zero and depend on the supcrlattice direction G. The 
constituent strain energy term Ai?cs( (J ) in Eq- (10) is 
specifically designed to reproduce these superlattice en- 
ergies, which are calculated directly from the LDA as 
follows. 

In the long-period limit pq — > oo the intcrfacial en- 
ergy becomes negligibly small (C(l/p)) in comparison 
with the elastic strain energy needed to deform the con- 
stituents to a common in-plane lattice constant a s . 55 ' 96 
Therefore, the formation energy per atom of Ar^B^ su- 
perlattice along G with composition x is given by the con- 
stituent strain energy AEcs{ x , G), defined as the equilib- 
rium (eq) value of the composition-weighted sum of the 



7 



Epitaxial parameters: Cu 



3 

V 

a 1 

s 
'3 

at 

a 

e 



0.5 
0.4 
0.3 
0.2 



W 0.1 



0.0 



fee Cu 


(110) ; 




Si) : 


(110)/ 
// (/ui; 


r-^° (ill) : 


(100) 
a 

i > 


(100) 

eq ^^-^ 

c , , , 



3.4 



3.6 



3.8 



4.0 



Substrate lattice constant a (A) 



FIG. 2. q(G) of fee Cu for principle directions as functions 
of the substrate lattice parameter a s . Directly calculated 
LDA values are represented by open symbols, lines show the 
fit using the expansion of 7(G) in Kubic harmonics. 



routinely used for semiconductor systems, 97 ' 100,101 
g's which do not depend on a s : 



9harm(G) = 1 — 



B 



Cu + A 7harm(G) 



gives 



(16) 



where 7harm(G) is a geometric function of the spherical 
angles formed by G: 



7harm(0, 0) = sin 2 (2#) + sin» sin 2 



(17) 



and Ki are the Kubic harmonics of angular momentum 
I. Figure 2 shows that the harmonic approximation man- 
ifestly breaks down for large epitaxial strains in met- 
als since there are several important qualitative differ- 
ences between the behavior in Fig. 2 and that predicted 
by the harmonic elasticity. First, q(a s ,G) strongly de- 
pends on the substrate lattice constant, while the har- 
monic q , harm(G') does not. Second, the harmonic expres- 
sion gives a definite order of q(G) as a function of the 
direction, i.e., cither (100) is the softest and then (111) 
must be the hardest, or vice versa. This order does not 
hold for large deformations. For instance, (201) becomes 
the softest direction for a s <C ao and (110) is the hardest 
for a s 3> a in Cu. Finally, g(100) exhibits a particularly 
dramatic softening for a s 3> ao, which has important con- 
sequences for the constituent strain energy and stability 
of superlattices along this direction. 96 

The above mentioned properties of qc u can be de- 
scribed by generalizing Eq. (17) for 7 to higher Kubic 
harmonics and strain-dependent expansion coefficients: 



j(a s ,G) = ^Tb l (a s )K l (G), 



(18) 



which has the property that in the harmonic limit 
(a s — > ao) all expansion coefficients with angular mo- 
menta higher than 4 tend to zero, reproducing 7h a rm from 
Eq. (17). Due to the cubic symmetry, only terms with 
I = 0,4, 6, 8, 10, 12, . . . enter in this expansion. Detailed 
discussion of the nonlinear epitaxial strain properties of 
elemental metals will be given in a separate publication. 96 

The constituent strain energy AE^? s (x,G) is calcu- 
lated numerically from Eq. (14) using the direct LDA 
values of AE cpl (a s ,G) for six principle directions. The 
obtained AE^ for these directions are shown in Fig. 3, 
illustrating several properties of the constituent strain 
which cannot be reproduced by the harmonic theory. 65 
First, the curves in Fig. 3 are skewed to different sides, 
while the harmonic Ais^g must be all skewed to the same 
side. Second, the calculated AE^? S cross for different di- 
rections, a property not allowed by the harmonic func- 
tional form. These crossings lead to (201) as the softest 
direction below x w 0.2, and (110) as the hardest for 
Au-rich superlattices, while the harmonic theory gives 
A£° q s (lll) as the highest and AE^(100) as the lowest 
constituent strain for all compositions of the studied no- 
ble metal alloys. The behavior of AE^? S for (100) is par- 
ticularly interesting, since the curves in Fig. 3 abruptly 
change slope around x w 0.15 and have very low values 
for x > j. As we show in Ref. 96, this is a manifestation 
of the low energy cost of deforming fee Cu into the body- 
centered tetragonal structure along the epitaxial Baincs 
path. Small constituent strain of (100) superlattices has 
profound influence on the predicted ground states of Cu- 
Au (see Sec. IV A 1). 

The constituent strain energy for arbitrary direction 
G is then obtained by interpolating between the prin- 
ciple directions using the following expansion in Kubic 
harmonics: 



AE CS (x,G) 



£>(aO#j(G). 
1=0 



(19) 



We have taken Z max = 10, which gives five composition- 
dependent fitting coefficients determined from a fit to 
the directly calculated values [Eq. (14)] for six princi- 
ple directions. The characteristic errors of this fit at 
the equiatomic composition are 1 — 2 meV/atom. Equa- 
tion (19) is then used in Eqs. (11)-(12). 



C. Constructing the Cluster Expansion 

Once we have a closed-form expresion for the equi- 
librium constituent strain energy AEcs(&) and a set 
{AiJ LDA (er)} of T = formation enthalpies, we deter- 
mine the unknown cluster interactions of Eq. (10) in the 
following two-step process: 
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FIG. 3. Equilibrium constituent strain energies for Cu-Au, 
Ni-Au and Cu-Ag. The constituent strain energy of Ag-Au is 
negligibly small and therefore not shown. 



First, the total energies of all structures from Table III 
are used in the fit to investigate the behavior of the root- 
mean-square (rms) error A rms of the fit, Eq. (13), as a 
function of the number of real-space pair and multibody 
interactions. Reciprocal space CE allows to add pair in- 
teractions systematically in the order of increasing in- 
tersite separation, up to any number of near-neighbor 
shells. The k-space smoothness criterion in Eq. (13) 
automatically selects optimally short-ranged interactions 
and chooses physically important pair interactions which 
are essential to produce a good fit to the directly calcu- 
lated LDA energies. The dependence of the rms error on 
the number of pair and multibody interactions is shown 
in Fig. 4. Figure 4(a) is obtained by fixing the number of 
multibody interactions, and varying the number of pair 
interactions. It shows that in all systems the cluster ex- 
pansion is well converged using 10 to 20 pair interactions. 
The convergence rate is fastest for Ag-Au and slowest for 
Ni-Au, which we attribute to increasing size mismatch 
going from Ag-Au to Ni-Au, with Cu-Ag and Cu-Au ex- 
hibiting intermediate convergence rates. 

Selection of important multibody interactions is more 
delicate. The number of pair interactions is fixed to a 
converged value (20 or more) , and a large set of 3- to 4- 
body figures is tested as to whether it improves the rms 
error of the overall fit. It is retained in the CE only if 
A rms decreases considerably. During the fitting process, 
we also monitor the overall stability of the CE, as mea- 
sured by a change in other multibody interactions upon 
the addition of a particular figure. Unstable behavior 
usually signals of linear dependencies in the chosen set 
of clusters and an ill-conditioned inverse problem, neces- 
sitating a different choice of {J/}. Figure 4(b) shows 
the convergence of the CE with respect to the number 
of multibody interactions, keeping N pa i rs equal to their 
converged values. An important thing to notice is that 
the multibody interactions produce a decrease in the rms 
error which is of the same magnitude as that due to the 
pair interactions. Furthermore, the effect of multibody 
interactions is largest in Ni-Au, and decreases in order of 
decreasing size mismatch, becoming negligible in Ag-Au. 

In the second step we test the stability of the fit and its 
predictive power. Using the trial set of figures obtained 
in the previous step, we exclude several structures which 
are fit rather well (e.g., Z2, (32, and LI2 in Ni-Au), and 
repeat the fit, obtaining new values of the effective cluster 
interactions. These values are used to predict the total 
energies of the structures excluded from the fit. If the 
change in AHce{o~) is not acceptable (more than few 
meV/atom), we return to the first step to search for a 
better set of interactions. The most severe test is to 
exclude structures with the poorest fit to their formation 
enthalpies, e.g., SQSU a and SQSUb in Ni-Au. If the 
predicted formation energy does not change significantly, 
the chosen set of figures is considered to be stable and 
predictive. The final cluster expansion is produced by 
using this set of figures and all structures from Tabic III. 

Figure 5 shows the calculated pair interactions as func- 
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FIG. 4. Root-mean-square errors A rms of the cluster expansions for Ag-Au, Cu-Ag, Cu-Au and Ni-Au as functions of the 
number of pair and multibody interactions. 



tion of the near-neighbor fee shell. There are several 
noteworthy trends in the four alloy systems: 

(i) Only in Ag-Au and Cu-Au are the nearest-neighbor 
pair interactions dominant: in Cu-Ag the 1-st and 3-rd 
neighbor pair interactions are of similar magnitude, while 
the 3-rd neighbor interaction dominates in Ni-Au. 

(ii) The dominant interactions have signs consistent 
with the observed phase diagrams: Ag-Au and Cu- 
Au have positive ( "antiferromagnetic" ) nearest-neighbor 
pair interactions J2, corresponding to the tendency to- 
wards complete miscibility and ordering at low tempera- 
tures. The behavior of Ni-Au, in spite of positive 1-st and 
2-nd neighbor pair interactions, is dominated by the "fer- 
romagnetic" 3-rd neighbor interaction L 2 (which causes 
phase separation at low temperatures). Both dominant 
1-st and 3-rd neighbor pair interactions in Cu-Ag are neg- 
ative, implying a miscibility gap. The constituent strain 
energy A£^,g is always positive and therefore increases 
the propensity for incoherent phase separation. 

(iii) Although the nearest-neighbor pair interaction is 
clearly dominant in Cu-Au, other pair interactions show a 
long-ranged oscillatory behavior extending over approx- 
imately 15 shells. As found in other systems, 65 ' 85 this 
is a direct consequence of the atomic relaxation caused 
by the constituent size mismatch between Cu and Au. 
The pair interactions are slowly decaying in Cu-Ag and 
Ni-Au, too. 

The calculated multibody interaction energies are 
shown in Figure 6. J\ is the point interaction, J3, K3, A3, 
are triplets and J4, K4, and L 4 are four-point clus- 
ters in increasing order of interatomic separation (see Lu 
et al. 54 for a full description of the clusters). Figure 6 
illustrates the importance of the multibody terms in our 
Hamiltonian. 



D. Finding the T = ground states and T > 
properties 

Having parametrized the configurational energies in 
terms of the mixed-space cluster expansion Eq. (10), 
we can use it with established statistical methods to 
predict various structural properties: T = ground 
states, order-disorder transition temperatures, configu- 
rational entropies, free energies, phase stabilities and 
atomic short-range order parameters. Due to the pres- 
ence of both reciprocal and real space terms in the Hamil- 
tonian (10), traditional techniques, e.g., the Cluster Vari- 
ation Method, are not readily applicable. Monte Carlo 
simulations must be used instead to calculate statistical 
properties at finite temperatures. The basic computa- 
tional algorithm is as follows. We adopt the Metropolis 
algorithm in the canonical ensemble (fixed composition). 
For each attempted spin flip, the change in the multiplct 
interaction energy is evaluated in the real space. To ob- 
tain the reciprocal space energy (constituent strain and 
pair interaction energies), the Fourier transform of the 
spin function S(Ri, a) is needed. It can be calculated ei- 
ther with the help of the Fast Fourier Transform (FFT) or 
evaluated directly taking advantage of the special method 
described in Ref. 87, which is much more economical: if 
the total number of sites in the simulation box is N, a full 
FFT has to be done only once after approximately every 
y/N accepted spin flips, which makes the whole compu- 
tational effort for this special method scale as N 1 - 5 . 

A simulation box of N = 4096 atoms (16 x 16 x 16) is 
used to calculate all thermodynamic properties presented 
in this paper. The transition temperatures are computed 
by cooling the system from high temperatures and moni- 
toring the discontinuities in the average energy and peaks 
in heat capacity. To eliminate possible hysteresis effects, 
the resulting low-temperature configurations are gradu- 
ally heated up past the transition point. The former 
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FIG. 5. Real space pair interactions for the studied noble 
metal alloy systems. 
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process provides the lower bound on the transition tem- 
perature, T\, while the latter gives the upper bound, T%. 
The heating and cooling rates are such that T\ and T 2 
differ by no more than 20K, an insignificant uncertainty 
compared to the inaccuracies of the LDA calculations and 
the fit errors of the cluster expansion. 1000 flips/site and 
a temperature decrease of 2% for each Monte Carlo step 
are usually sufficient, although in a few cases the results 
are checked using 2000 flips/site and 0.5% temperature 
change. 

Zero temperature ground states are found by cooling 
the system to T = and checking whether the energy 
of the final configuration lies on the convex hull. This 
process is repeated for several random number seeds 
and starting temperatures, always yielding configura- 
tions with similar (usually identical) energies. We ex- 
plore many equally spaced compositions with an interval 
Ax = 0.05. The number of possible configurations for 
each x is N coni = ( xj v)t(iv(i-x))! • 

Configurational entropy of the disordered alloys at fi- 
nite T is computed from the energy vs. temperature 
curves obtained by cooling the system from very high 
("T = oo") temperatures. The following thermodynamic 
formula gives the configurational entropy at temperature 
T: 

AS coni (T) = A5 idoal + E(T)/T - k B / E{(3') dp 1 , 

Jo 

(20) 

where (3 = l/fc^T and ASideai = fce^loga; + (1 — 
x) log(l — x)] is the configurational entropy of an ideal 
solid solution. 



IV. RESULTS 
A. T — Ground States 

1. Ground states of Cu-Au 

Figure 7 shows the calculated T = ground state lines 
of Cu-Au and Ag-Au which were obtained from simulated 
annealing quenches of a 16 x 16 x 16 system. In Cu-Au, 
we find the Ll 2 (CU3AU) and LIq (CuAu) structures as 
the stable ground states of Cu-rich alloys, in agreement 
with the existing phase diagram data. 1-4 These data also 
list LI2 as the stable low-temperature phase of CuAu3. 
However, we find new, previously unsuspected ground 
states of Au-rich compounds, all belonging to the fam- 
ily of (001) superlattices. At x — | we find a stable (32 
(CuAu2) phase (prototype MoSi2), which is a CuiAu2 su- 
perlattice along (001). At x = §, our cluster expansion 
predicts that a complex CU1AU4CU1 AU4CU1AU2CU1AU2 
(001) superlattice falls on the convex hull, although its 
energy is less than 2 meV below the tieline connecting 
(32 (CuAu 2 ) and Au. Furthermore, even the directly 
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Ag Composition x Au 

FIG. 7. T = K ground state lines for Cu-Au and Ag-Au 
obtained from simulated annealing calculations. LI2 CuAu3 
is not only above the ground state line, but also has a higher 
formation enthalpy than other structures at the same com- 
position, e.g. LDA calculations place the formation enthalpy 
of Z3 below that of LI2. Plots for Cu-Ag and Ni-Au are not 
shown since these systems phase separate at T = K. 

calculated LDA enthalpy of formation of Z3 (which is 
a CU1AU3 (001) superlattice) is considerably lower than 
that of LI2 CuAu3. 

We carefully checked whether the predicted new LDA 
ground states for Au-rich Cu-Au alloys artifacts of some 
approximation in our LDA calculations or the fit error of 
the cluster expansion. The latter possibility was quickly 
dismissed, since the directly calculated LDA enthalpies 
of formation for LIq, (32, Ll 2 and Z3 agreed with the 
values derived from the cluster expansion to better than 
2 mcV/atom (see Table III), while the new (100) SL 
ground state is 14 meV/atom below Ll 2 . To address 
the former possibility, we performed careful convergence 
tests for LIq, (32, Ll 2 and Z3 with respect to the plane 
wave cutoff and number of k points in the first Bril- 
louin zone. The cutoff was increased from i?_ftT max = 9 
to RK max = 11 and the density of the Brillouin zone 
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mesh was doubled from 8 x 8 x 8 to 16 x 16 x 16, an 
eightfold increase in the total number of k points. These 
tests showed that the formation enthalpies of Llo, (32 
and Ll 2 were converged to within 1 meV/atom with re- 
spect to the size of the basis set and the number of k 
points. Further, we checked how the choice of muffin- 
tin radii affected AH. Varying 7?mt(Au) between 2.3ao 
and 2.5ao changed the formation enthalpies by at most 
2 meV/atom and did not shift the relative stabilities of 
phases. Finally, we repeated these calculations using 
the Perdew-Zunger 90 parametrization of the Ceperley- 
Aldcr 91 LDA functional, as well as the generalized gra- 
dient approximation (GGA) of Perdew and Wang, 92 and 
found insignificant (about 2 meV/atom) changes in the 
formation enthalpies. Inclusion of the spin-orbit inter- 
action in the second variation procedure 104 changed the 
formation enthalpy of LIq (CuAu) by only 3.7 meV/atom 
(from —48.2 to —51.9), indicating that it is not im- 
portant for the energetics of Cu-Au. This conclusion 
is in line with the findings of Ref. 105 that the spin- 
orbit interaction influences the band structure but has 
little effect on equilibrium lattice properties. Therefore, 
we conclude that state-of-the-art first-principles density 
functional calculations do not predict Ll 2 to be a stable 
T = ground state of CuAu^. It is possible that van 
der Waals interactions, omitted by the LDA and impor- 
tant for large, polarizable atoms such as Au, can affect 
the formation energies and hence the ground states of 
Cu-Au. 

We next analyze the possibility that the correct T = 
ground state around x = | is not Ll 2 as has 
been assumed in the literature before. Although most 
compilations 1-4 of binary alloy phase diagrams give Ll 2 
as the stable structure of CuAu3, the experimental 
evidence 7,8,10 seems inconclusive because of the difficul- 
ties in obtaining equilibrated long-range ordered samples. 
X-ray studies 8 have found superlattice peaks consistent 
with the cubic Ll 2 structure, but only very broad low- 
order reflections have been observed. These superlattice 
lines could not be sharpened by any heat treatment. 8 It 
is not clear to us if the X-ray reflections can be rein- 
dexed according to some other non-Ll2 phase. It is also 
possible that at elevated (T 500 K) temperatures Ll 2 
is stabilized by the entropy (configurational and vibra- 
tional), while another transformation to the low-energy 
structure should occur but is kinetically inhibited below 
500 K. The biggest experimental obstacles to verifying 
our predictions seem to be low diffusion rates below the 
ordering temperature of CuAu3, T c 500 K. 

Next we discuss the experimental signatures of the new 
LDA ground state structures. MoSi2-type /32 CuAu2 has 
a superlattice reflection at (§00), but CuAu 3 (100) su- 
perlattice has reflections at (100) and (|00). These re- 
flections also manifest themselves in the predicted atomic 
short-range order of the disordered alloys (for details see 
Ref. 103). 



2. Ground states of Ag-Au, Cu-Ag and Ni-Au 

The ground state line of Ag-Au is shown in Figure 7(b), 
exhibiting Ll 2 (Ag 3 Au), Ll a (AgAu) and Ll 2 (AgAu 3 ) 
stable low-temperature phases. Experimentally, these al- 
loys are known to be completely miscible, 2-4 and there 
are several indications 69 that they would order below 200 
K if not for the very low diffusion rates. Theoretical 
transition temperatures and short-range order patterns, 
as well as a complete discussion are given by Lu and 
Zunger. 54 

The calculated ground states of Cu-Ag and Ni-Au are 
found to be phase separation, in agreement with the ex- 
perimental enthalpy data. 2 Neither alloy has a single or- 
dered or disordered structure with negative enthalpy of 
formation and therefore there are no stable T — ground 
states except the phase-separated alloy. 

B. Mixing enthalpies 

It is interesting to compare the calculated mixing en- 
thalpies of disordered Cu-Au alloys with the available 
theoretical and experimental data. Table IV summarizes 
the values of AH m i x (x,T) for the completely random 
(T = 00), short-range ordered (T = 800 K) and com- 
pletely ordered (T = K) Cu-Au alloys at compositions 
x = 3' \ an(1 f ■ Several important points are apparent 
from this table: 

(i) Studies 50,48,62 which have completely neglected 
atomic relaxations predict a substantially positive en- 
thalpy of formation for the completely random alloy. In 
our calculations, relaxations in the random alloy reduce 
AH m i x (T = 00) by a large amount, bringing it down to 
essentially zero. 

(ii) Comparison of the present results for the T = 00 
random alloys with those of Wei et a/. 51 shows the in- 
fluence of the number of structures included in the clus- 
ter expansion. Since Wei et al. used the same FLAPW 
method 88 , but included a set of only five high-symmetry 
ordered structures [Al (Cu), Ll 2 (Cu 3 Au), Ll (CuAu), 
LI2 (CuAu 3 ) and Al (Au)], the atomic relaxation effects 
were included incompletely. Indeed, their treatment gives 
much larger mixing enthalpies of the random Cu-Au al- 
loys than the present work employing approximately 30 
low-symmetry structures with large relaxations. There- 
fore we conclude that the Connolly- Williams set of five 
ordered structures cannot correctly capture the large de- 
crease of the mixing enthalpy of random Cu-Au alloys 
caused by the atomic relaxations. 

(iii) The good agreement between the relaxed (this 
study) and "unrelaxed" (Wei et al. 51 ) values of Ai? m ; x 
at T = 800 K suggests that the short-range order in Cu- 
Au tends to decrease the role of the atomic relaxations. 
This effect can be qualitatively explained on the basis of 
the ordering tendency towards high-symmetry structures 
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TABLE IV. Calculated mixing enthalpies of disordered Cu i _ I Au I alloys compared with the values obtained by other studies 
and experimental measurements (in meV/atom). FLAPW is the full-potential linearized augmented plane wave method, LMTO 
- linearized mumn-tin-orbitals method, KKR - Korringa-Kohn-Rostoker multiple scattering method, ASA - atomic-sphere 
approximation, CPA - coherent potential approximation, CWM - Connolly- Williams cluster expansion, MSCE - mixed-space 



cluster expansion 


used in this 


study, "Rel." - 


incorporating atomic relaxations, and 


"Unrel." 


- neglecting atomic relaxations. 


Composition 


Expt. f 


This 


Wei 


Amador Terakura 


Ruban 


Weinberger 








et at. 


et aL. 


et at. 


et al. 


et aL. 






FLAPW 


FLAPW 


T ft/TT^O A Q A 


Ao VV 


LMTO-ASA 


KKR-ASA 






MSCE 


CWM 


W 1V1 


O VV 1V1 


CPA 


CPA 






(Rel.) 


(Rel.) 


(Unrol.) 


(Unrol.) 


(Unrcl.) 


(Unrcl.) 


AH ■ (T = oo) 
















Cuo,75 Alio. 25 




-\-z.o 


+46.3 


+59 


+26.9 


+54.6 


— 27 


Cuo.5oAuq.50 




+1.6 


+38.0 


+61 


+30.4 


+44.3 


-57 


Cu 25 Au 0.75 




+5.4 


+ 18.6 


+39 


+20.4 


+ 19.8 


-31 


A« mi x(T = 800 K) 
















Cuo.75Auo.25 


-46 g 


-17.3 




-6 








Cuo.5oAuq.5o 


-53 g 


-19.3 


-16.9 


-5 








Cuo.25Auo.75 


-31 s 


-1.2 


-2.6 


+8 








AH mix (T = K) 
















Ll 2 Cu 3 Au 


-74 


-37.3 


-36.0 




-65.0 


-60.7 


-54 


L1q CuAu 


-91 


-48.2 


-62.9 




-69.7 


-83.4 


-76 


LI2 CuAu3 


-59 


-17.3 


-26.4 




-34.0 


-56.1 


-47 



a Ref. 51 using the Connolly-Williams structures (relaxation of Ll only). 

b Ref. 50. 

c Ref. 48. 

d Ref. 62. 

c Ref. 58. 

f Ref. 2. 

s Values obtained at T = 720K. 



which have little or no relaxation energy (LI 2 and L lo- 
in Cu-rich alloys). 

(iv) The mixing enthalpies of the random alloy calcu- 
lated by Weinberger et al. 58 using the coherent-potential 
approximation (CPA) differ strongly not only from those 
obtained using the cluster expansion methods, 51 ' 50,48 but 
also from the numbers given in the CPA work of Ruban, 
Abrikosov, and Skriver. 62 Since the CPA of Weinberger et 
al. 58 neglects the (a) atomic relaxation, (b) charge trans- 
fer and (c) short-range order, which all lower the forma- 
tion energies, the negative values obtained by Weinberger 
et al. 58 are very puzzling. 

(v) There are significant discrepancies between the best 
calculated and experimentally measured 15,1,2 values of 
AiJ mix at both T = K and T — 800 K. At present 
these discrepancies are hard to explain since the available 
general potential LDA calculations 51,52,57 of AH(L1 2 ) 
and AH(L1 ) agree with each other reasonably well. On 
the other hand, formation energies in Cu-Au are numer- 
ically very small and present a severe test for any first- 
principles model of electronic exchange-correlation. It is 
noteworthy that several less accurate first-principles cal- 
culations, using the atomic-sphere approximation (ASA), 
have achieved better agreement with the experimental 
enthalpies of formation than the state of the art gen- 
eral potential techniques. We consider this to be fortu- 
itous. In all cases, LDA calculations correctly predict the 
relative magnitudes of AH for Ll 2 and Ll , as well as 
reproduce measured asymmetry in formation enthalpies 
towards more negative values of AiJ m ; x for Cu-rich al- 



loys. 

C. Order-disorder transition temperatures 

Order-disorder transitions have been investigated at 
compositions (x = j, i, | and |) using the Monte Carlo 
simulation technique described in Sec. HID. The result- 
ing transition temperatures, T c , are given in Table V. 
All transitions are found to be first order, involving dis- 
continuities in the energy and correlation functions. At 
x = j we find a transition from the disordered state to 
long-range ordered Ll 2 Cu 3 Au at T c = 530 K, which is 
only 130 K lower than the experimentally observed tran- 
sition temperature. For the equiatomic alloy at x = | 
the calculated and experimental transition temperatures 
agree to a few degrees Kelvin. However, we do not find 
the CuAu II phase which exists in a narrow tempera- 
ture range between 658 K and 683 K. This phase is sta- 
bilized by the free energy differences between LIq and 
long-period superstructures of Ll which are as small as 
1 meV/atom 56 and therefore beyond the accuracy of self- 
consistent LDA calculations. 

For x — I we obtain a sequence of transformations, 
the first one occuring at T = 750 K from the disordered 
Al phase to a coherent two-phase mixture of 02 and Al. 
Then a subsequent transition at T — 635 K takes CuAu 3 
into the long-range ordered (100) superlattice which is 
predicted to be the stable T = ground state at that 
composition (sec Sec. IVA1). The calculated transition 
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TABLE VI. The experimentally measured 2 entropy of formation ASl°l m , the calculated configurationl entropy AS^i an d 
the derived non-configurational entropy of formation, AS n °7_ coll f ■• All values are given in units of fee/atom. 



System 


X 


T(K) 


A oform 
'-"tat 


AS'ideal 


A ocalc 
^■"conf 


A oform 

non — conf 

AStl m - 


Cu-Au 


0.5 


800 


0.73 


0.69 


0.57 


0.16 


Ag-Au 


0.5 


800 


0.52 


0.69 


0.62 


-0.10 


Cu-Ag 


0.141 


1052 


0.77 


0.41 


0.40 a 


0.37 


Ni-Au 


0.5 


1100 


1.04 


0.69 


0.56 


0.48 



a This value was obtained at T = 1136 K, since a coherent phase separation starts at lower temperatures. 



TABLE V. Calculated order-disorder transition tempera- 
tures (in K) for Cu-Au. Al denotes the configurationally dis- 
ordered fee phase, and n/a means that the transition has not 
been observed (either experimentally or in the Monte Carlo 



simulation) . 


Compo- 


Tran- 


Expt. 


This 


sition 


sition 




study 


Cuo.75 Auo.25 


Al -f Ll 2 


663 


530 


Cuo.50 Auo.50 


Al -» Llo 


683/658 a 


660 


Cuo.33 Auo.66 


Al -> (32 


n/a 


735 


Cuo.25 Auq.75 


Al -f Ll 2 


« 500 


n/a 




Al (32 + Al 


n/a 


750 




(32 + Al^ (100)SL 


n/a 


680 



a CuAu undergoes a transition to CuAu-II at 683 K, subse- 
quently transforming into Llo CuAu-I at 658 K. 



at x = | goes straight into the (32 phase at T = 735 K. 
Therefore, a two-phase (32 + Al field is predicted to exist 
at temperatures somewhere between 635 K and 730 K 
and around x = |. These predictions reflect the LDA. 
As stated in Sec. IV A 1, corrections to the LDA might 
be significant. 

D. Non-configurational entropy 

The effect of the non-configurational entropy (elec- 
tronic, vibrational, etc.) on the alloy phase stability 
has recently attracted considerable interest. 106-116 For 
instance, it has been suggested 108-115 that there are 
large differences in the vibrational entropies of order- 
ing Adored — •S'disord' wmcn should manifest themselves 
in shifts of the order-disorder transition temperatures. 
There is another important class of thermodynamic prop- 
erties where the vibrational entropy may play a role, and 
which has often been overlooked. Namely, it is the en- 
tropy of formation with respect to the pure constituents, 
defined in analogy with AH in Eq. (4): 

AStt m (M- x B x ,T) = S(A^ X B X ,T) 

— (1 - x)S(A, T) - xS(B, T), (21) 

where S(A,T) is the total entropy of the pure con- 
stituent A at temperature T. It is often assumed that the 
configurational entropy is the dominant contribution to 



ASl° T t m (Ai- x B X} T) because all other contributions can- 
cel out in Eq. (21). The non-configurational entropy of 
formation, 

A^ r n m -conf(Ai-,B x ,T) = AS t f °r(Ai-,B x ,T) 

-A5 conf (A 1 _ ;c B x ,T), (22) 

contributes to such important quantities as mutual solu- 
bility limits and miscibility gap temperatures. 

Noble metal alloys are excellent cases to test the val- 
ues of AS£°™ 

conf since accurate experimental data on 
the entropies of formation, AS*£™, are available, and 
the configurational entropy AS COIL { can be calculated 
accurately using the thermodynamic integration tech- 
nique described in Sec. HID. Table VI gives the mea- 
sured entropies of formation for disordered solid solutions 
Ai-zB-e, ASl°l m (x, T), the maximum attainable config- 
urational entropy AS'ideal, as well as the theoretically cal- 
culated configurational entropy A5^ c f , and the derived 
value for the non-configurational entropy of formation, 
conf . It shows that the size-mismatched noble 
metal systems have large amounts of AS^"™ conf in the 
disordered solid solution. Since it is unlikely that these 
values of AiS^Jj^l f are of electronic or magnetic origin, 
we suggest that the excess entropy in the disordered solid 
solutions of Ni-Au, Cu-Ag and Cu-Au is vibrational. It 
is possible that the atomic relaxations lead to a softening 
of lattice vibrations, although the physical mechanism of 
this softening is unclear at present. 

Sanchez et al. 49 in their study of the Cu-Ag system 
noted that even a very crude model of the vibrational 
entropy markedly improved the agreement with the ex- 
perimental solubility data. In the case of Ni-Au, which 
exhibits the largest A5^™ conf , it is possible to reconcile 
the experimentally measured and theoretically calculated 
miscibility gap temperatures only by taking into account 
the non-configurational entropy of formation. 117 

The fact that Cu-Au also has a positive AS^^_ con{ has 
little qualitative effect on the phase diagram since Cu 
and Au are completely miscible from total energy and 
configurational entropy considerations alone. Ag-Au is 
calculated to have a negative A5^°™ conf , but its value is 
close to the experimental uncertainty in the measurement 
of AS. 
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FIG. 8. SQS bondlengths for Cu-Au and Ni-Au. 

E. Bond lengths in random alloys 

Since recent experimental measurements of the 
composition-dependence of interatomic bond lengths in 
Cu-Au 24 and Ni-Au 23 have found several unusual fea- 
tures, it is interesting to address these trends from first- 
principles LDA calculations. In the present work we 
model the atomic positions in the random alloys using 
special quasirandom structures 118 (SQS). These periodic 
structures are designed to reproduced the pair and multi- 
body correlation functions of the perfectly disordered 
configuration as closely as possible. It has been shown 118 
that even small unit cell SQS's can give rather accurate 
representation of the properties of random alloys. We 
have performed LDA calculations for 8 atom/cell SQS's 
at x = \ (SQSU a ), x = i (SQS8 a ,SQS8 b ) and x = § 
(SQSlAb). The atomic positions and cell coordinates 



have been fully relaxed to minimize the total energy. The 
results for Cu-Au and Ni-Au interatomic bond lengths 
are shown in Fig. 8. The main features are: 

(i) In spite of the different phase diagram properties 
(Ni-Au phase separates and Cu-Au orders at T = K), 
the calculated behavior of bond lengths is very similar, 
which we attribute to the similar size mismatch in both 
systems (12% in Cu-Au and 15% in Ni-Au). 

(ii) Our calculations give three distinct bond 
lengths at all compositions, which is also observed 
experimentally 23 ^ 24 Probably the most interesting fea- 
ture in Fig. 8 is the crossing of Rbb(x) and Rab{x) 
curves at x = | in both systems. The measurements 
for Cu-Au 24 and Ni-Au 23 indicate that this may indeed 
be correct, since the deduced values around this compo- 
sition are very close and have large error bars. 

(iii) Another important feature, observed experimen- 
tally and reproduced by our SQS results, is that A — A 
bonds change much more as x varies from to 1 than 
B — B bonds when x varies from 1 to 0, suggesting that 
the compressed bonds become increasingly stiff and the 
expanded bonds weaken. This behavior can be explained 
by the asymmetry in the interatomic potential curves, 
which are rapidly hardening upon compression and soft- 
ening upon expansion. However, our results for Raa at 
x = | and fisB at i = | are obtained from an aver- 
age of only 4 minority bonds in the SQS1A structures, 
and perhaps are not representative of a wider statistical 
sample. 

(iv) It is interesting to note that the predicted bond 
lengths between unlike atoms Rab do not follow the lin- 
ear relation Rab = Raa + x(Rbb — Raa)- 



V. SUMMARY 

We have showed that accurate first-principles stud- 
ies of alloys with large size mismatches are now feasible 
using the mixed-space cluster expansion method. This 
method has been applied to noble metal alloys where 
vast amounts of experimental data and many theoretical 
studies are available. 

(i) The mixed-space cluster expansion has been gen- 
eralized to include the effects of nonlinear strain on the 
formation energies of long-period superlattices. We find 
that the elastic energy, required to lattice-match Cu and 
Ni to (100) surfaces of Au and Ag, is anomalously low, 
leading to a very low constituent strain energy of (100) 
superlattices. This effect is partly responsible for the sta- 
bilization of new LDA ground states of Au-rich Cu-Au 
alloys. 

(ii) In Au-rich Cu-Au, we predict new T = K 
ground states. Our LDA results place LI2 (CuAus), 
previously thought of as the stable T = state of 
CuAu3, higher in energy than a family of superlat- 
tices along (100) direction. In particular, MoPt2-type 
CuAu2 [CU1AU2 superlattice along (100)] and a compli- 
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cated CuiAu4CuiAu4CuiAu 2 CuiAu 2 (100) superlattice 
are found to be the LDA ground states. 

(hi) There are significant discrepancies (up to 50%) be- 
tween the experimentally measured and calculated LDA 
mixing enthalpies for Cu-Au alloys. This is surprising 
since the experimental mixing enthalpies of Ni-Au and 
Ag-Au are reproduced very well. 54 ' 117 

(iv) The calculated order-disorder transition temper- 
atures are in an excellent agreement with experiment. 
For instance, T c calc (x = ±) = 530 K and T c calc (x = 
\) = 660 K, compared with T^{x = \) = 663 K and 
T c ex P* (x = ±) = 683/658 K. 

(v) From the experimentally measured entropies of 
formation ASl°l m and the calculated configurational 
entropies AS^^, we obtain large non-conhgurational 
(probably vibrational) entropies of formation in the size- 
mismatched systems, Afc conf = AS t f ° t rm - AS C C ^ { . 
These entropies allow one to reconcile the experimental 
miscibility gap temperature and formation enthalpies of 
Ni-Au with the theoretical LDA values. 117 

(vi) Bond length distributions in Ni-Au and Cu- 
Au have been studied via supercell calculations em- 
ploying the special quasirandom structure technique. 
The important qualitative features of recent EXAFS 
measurements 23 ' 24 are correctly reproduced: existence of 
distinct A — A, B — B and A — B bond lengths at all 
compositions, possible crossing of Raa{x) and Rab(x) 
around x = | (where x is the composition of the larger 
constituent), softening of the shorter bond as i — ► 1, 
and deviations of the bond length Rab{x) between un- 
like atoms from the linear Vergard's law. 
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